// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtRevolveFilter.C
// ************************************************************************* //

#include <avtRevolveFilter.h>

#include <math.h>
#include <float.h> // for FLT_MAX

#include <vtkCell.h>
#include <vtkCellData.h>
#include <vtkIdList.h>
#include <vtkMatrix4x4.h>
#include <vtkPointData.h>
#include <vtkTemplateAliasMacro.h>
#include <vtkUnstructuredGrid.h>

#include <avtExtents.h>

#include <BadVectorException.h>
#include <InvalidCellTypeException.h>
#include <InvalidDimensionsException.h>
#include <vtkVisItUtility.h>

static void GetRotationMatrix(double angle, double axis[3], vtkMatrix4x4 *mat, 
                              avtMeshCoordType mt);


// ****************************************************************************
//  Method: avtRevolveFilter constructor
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Wed Dec 11 11:31:52 PDT 2002
//
// ****************************************************************************

avtRevolveFilter::avtRevolveFilter()
{
}


// ****************************************************************************
//  Method: avtRevolveFilter destructor
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Wed Dec 11 11:31:52 PDT 2002
//
//  Modifications:
//
// ****************************************************************************

avtRevolveFilter::~avtRevolveFilter()
{
}


// ****************************************************************************
//  Method:  avtRevolveFilter::Create
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Wed Dec 11 11:31:52 PDT 2002
//
// ****************************************************************************

avtFilter *
avtRevolveFilter::Create()
{
    return new avtRevolveFilter();
}


// ****************************************************************************
//  Method:      avtRevolveFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Wed Dec 11 11:31:52 PDT 2002
//
//  Modifications:
//    Kathleen Bonnell, Wed May 21 11:10:58 PDT 2003   
//    Test for bad Rotation Axis (0, 0, 0), throw exception when encountered.
//
// ****************************************************************************

void
avtRevolveFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const RevolveAttributes*)a;
    const double *axis = atts.GetAxis();
    if (axis[0] == 0. && axis[1] == 0. && axis[2] == 0.)
    {
        EXCEPTION1(BadVectorException, "Axis of Revolution");
        return;
    }
}


// ****************************************************************************
//  Method: avtRevolveFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtRevolveFilter with the given
//      parameters would result in an equivalent avtRevolveFilter.
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Wed Dec 11 11:31:52 PDT 2002
//
// ****************************************************************************

bool
avtRevolveFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(RevolveAttributes*)a);
}


// ****************************************************************************
//  Method: avtRevolveFilter::GetMeshType
//
//  Purpose:
//      Determine the mesh type.  Uses the data attributes unless overridden
//      by the operator attributes.
//
//  Programmer: Hank Childs
//  Creation:   March 18, 2007
//
// ****************************************************************************

avtMeshCoordType
avtRevolveFilter::GetMeshType(void)
{
    switch (atts.GetMeshType())
    {
      case RevolveAttributes::XY:
        return AVT_XY;
      case RevolveAttributes::RZ:
        return AVT_RZ;
      case RevolveAttributes::ZR:
        return AVT_ZR;
      default:
        break;   // fall thru to below
    }

    // Auto ... "VisIt Decides"
    return GetInput()->GetInfo().GetAttributes().GetMeshCoordType();
}


// ****************************************************************************
//  Method: avtRevolveFilter::GetAxis
//
//  Purpose:
//      Determines the axis.  
//
//  Programmer: Hank Childs
//  Creation:   March 18, 2007
//
//  Modifications:
//    Brad Whitlock, Mon Mar 30 10:00:45 PDT 2009
//    Added a missing return. I also added a default case.
//
// ****************************************************************************

void
avtRevolveFilter::GetAxis(avtMeshCoordType mt, double *axis)
{
    if (!atts.GetAutoAxis())
    {
        axis[0] = atts.GetAxis()[0];
        axis[1] = atts.GetAxis()[1];
        axis[2] = atts.GetAxis()[2];
        return;
    }
    switch (mt)
    {
    default:
    case AVT_XY:
        axis[0] = 1.;
        axis[1] = 0.;
        axis[2] = 0.;
        break;
    case AVT_RZ:
        axis[0] = 0.;
        axis[1] = 0.;
        axis[2] = 1.;
        break;
    case AVT_ZR:
        axis[0] = 0.;
        axis[1] = 0.;
        axis[2] = 1.;
        break;
    }
}

//
// Templated methods for revolving points and vectors
//

template <class T> void
RevolvePoints(int niter, int nsteps, vtkIdType npts, double stop_angle, 
    double start_angle, double axis[3], vtkMatrix4x4 *mat, avtMeshCoordType mt,
    vtkDataSet *in_ds, T *ptr)
{
    for (int i = 0 ; i < niter ; ++i)
    {
        double angle = ((stop_angle-start_angle)*i)/(nsteps-1) + start_angle;
        GetRotationMatrix(angle, axis, mat, mt);
        for (vtkIdType j = 0 ; j < npts ; ++j)
        {
            double pt[4];
            in_ds->GetPoint(j, pt);
            pt[3] = 1.;
            double outpt[4];
            mat->MultiplyPoint(pt, outpt);
            ptr[0] = (T)outpt[0];
            ptr[1] = (T)outpt[1];
            ptr[2] = (T)outpt[2];
            ptr+=3;
        }
    }
}

void
RevolvePoints(vtkDataArray *pts, int niter, int nsteps, vtkIdType npts, 
    double stop_angle, double start_angle, double axis[3], vtkMatrix4x4 *mat, 
    avtMeshCoordType mt, vtkDataSet *in_ds)
{
    switch(pts->GetDataType())
    {
        vtkTemplateAliasMacro(
            RevolvePoints(niter, nsteps, npts, stop_angle, start_angle, 
                          axis, mat, mt, in_ds,
                          static_cast<VTK_TT *>(pts->GetVoidPointer(0))));
    }
}

template<class T> void
RevolveNodeVectors(int istop, vtkIdType jstop, int nsteps, 
    double stop_angle, double start_angle, double axis[3], 
    vtkMatrix4x4 *mat, avtMeshCoordType mt, T* vecptr)
{
    for (int i = 0 ; i < istop ; ++i)
    {
        double angle = ((stop_angle-start_angle)*i)/(nsteps-1) + start_angle;
        GetRotationMatrix(angle, axis, mat, mt);

        for (vtkIdType j = 0 ; j < jstop ; ++j)
        {
            double in_pt[4];
            in_pt[0] = (double)vecptr[0];
            in_pt[1] = (double)vecptr[1];
            in_pt[2] = (double)vecptr[2];
            in_pt[3] = 1.;
            double out_pt[4];
            mat->MultiplyPoint(in_pt, out_pt);
            vecptr[0] = (T)out_pt[0];
            vecptr[1] = (T)out_pt[1];
            vecptr[2] = (T)out_pt[2];
            vecptr += 3;
        }
    }
}

template<class T> void
RevolveCellVectors(int istop, vtkIdType jstop, int nsteps, 
    double stop_angle, double start_angle, double axis[3], 
    vtkMatrix4x4 *mat, avtMeshCoordType mt, T* vecptr)
{
    for (int i = 0 ; i < istop ; ++i)
    {
        double angle = ((stop_angle-start_angle)*(i+0.5))/(nsteps-1) + 
                         start_angle;
        GetRotationMatrix(angle, axis, mat, mt);

        for (vtkIdType j = 0 ; j < jstop ; ++j)
        {
            double in_pt[4];
            in_pt[0] = (double)vecptr[0];
            in_pt[1] = (double)vecptr[1];
            in_pt[2] = (double)vecptr[2];
            in_pt[3] = 1.;
            double out_pt[4];
            mat->MultiplyPoint(in_pt, out_pt);
            vecptr[0] = (T)out_pt[0];
            vecptr[1] = (T)out_pt[1];
            vecptr[2] = (T)out_pt[2];
            vecptr += 3;
        }
    }
}

void
RevolveVectors(vtkDataArray *vectors, bool cellVectors, int istop, 
    vtkIdType jstop, int nsteps, double stop_angle, double start_angle, 
    double axis[3], vtkMatrix4x4 *mat, avtMeshCoordType mt)
{
    if (!cellVectors)
    {
      switch(vectors->GetDataType())
      {
        vtkTemplateAliasMacro(
            RevolveNodeVectors(istop, jstop, nsteps,
                           stop_angle, start_angle, axis, mat, mt, 
                           static_cast<VTK_TT *>(vectors->GetVoidPointer(0))));
      }
    }
    else 
    {
      switch(vectors->GetDataType())
      {
        vtkTemplateAliasMacro(
            RevolveCellVectors(istop, jstop, nsteps,
                           stop_angle, start_angle, axis, mat, mt, 
                           static_cast<VTK_TT *>(vectors->GetVoidPointer(0))));
      }
    }
}


// ****************************************************************************
//  Method: avtRevolveFilter::ExecuteData
//
//  Purpose:
//      Sends the specified input and output through the Revolve filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Wed Dec 11 11:31:52 PDT 2002
//
//  Modifications:
//
//    Hank Childs, Thu Mar  8 11:08:53 PST 2007
//    Initialize the transparency actor with a good default.  If there should
//    be a scale factor, that will be determined later.
//
//    Hank Childs, Sun Mar 18 10:12:57 PDT 2007
//    Add support for RZ and ZR meshes.
//
//    Brad Whitlock, Fri Oct 12 14:43:23 PST 2007
//    Added support for revolving lines.
//
//    Eric Brugger, Tue Jul 10 11:37:19 PDT 2012
//    Added support for vector data.
//
//    Kathleen Biagas, Wed Aug  8 15:05:26 PDT 2012
//    Support double-precision by use of templated methods for points and 
//    vectors.
//
//    Eric Brugger, Thu Aug 14 09:05:33 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
//    Alister Maguire, Fri Jun 19 12:06:51 PDT 2020
//    Added ability to revolve points (VTK_VERTEX).
//
// ****************************************************************************

avtDataRepresentation *
avtRevolveFilter::ExecuteData(avtDataRepresentation *in_dr)
{
    //
    // Get the VTK data set.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();

    // Get information about how to revolve.
    //
    int    nsteps      = atts.GetSteps();
    double start_angle = atts.GetStartAngle();
    double stop_angle  = atts.GetStopAngle();
    if (start_angle == stop_angle)
    {
        return in_dr;
    }
    if (start_angle > stop_angle)
    {
        start_angle = atts.GetStopAngle();
        stop_angle  = atts.GetStartAngle();
    }
    avtMeshCoordType mt = GetMeshType();
    double axis[3];
    GetAxis(mt, axis);

    //
    // Set up our VTK structures.
    //
    vtkMatrix4x4 *mat = vtkMatrix4x4::New();
    vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
    vtkPoints *pts = vtkVisItUtility::NewPoints(in_ds);
    vtkIdType npts = in_ds->GetNumberOfPoints();
    vtkIdType ncells = in_ds->GetNumberOfCells();
    vtkIdType n_out_pts = npts*nsteps;
    pts->SetNumberOfPoints(n_out_pts);
    ugrid->SetPoints(pts);

    //
    // Create the points for each timestep.
    //
    int niter = (fabs(stop_angle-start_angle-360.) < 0.001 ? nsteps-1 : nsteps);
    RevolvePoints(pts->GetData(), niter, nsteps, npts, stop_angle, start_angle, 
                  axis, mat, mt, in_ds);

    //
    // Now set up the connectivity.  The output will consist of revolved
    // quads (-> hexes), revolved triangles (-> wedges) and revolved points
    // (vertexes). No special care is given to the case where an edge of a
    // cell lies directly on the axis of revolution (ie: you get a degenerate
    // hex, not a wedge).
    //
    vtkIdType n_out_cells = ncells*(nsteps-1);
    ugrid->Allocate(8*n_out_cells);
    bool overlap_ends = (fabs(stop_angle-start_angle-360.) < 0.001);
    for (vtkIdType i = 0 ; i < ncells ; ++i)
    {
         vtkCell *cell = in_ds->GetCell(i);
         int c = cell->GetCellType();
         if (c != VTK_QUAD && c != VTK_TRIANGLE && c != VTK_PIXEL &&
             c != VTK_LINE && c != VTK_VERTEX)
         {
             EXCEPTION1(InvalidCellTypeException, "anything but lines, quads,"
                                                   " points, and tris.");
         }
         vtkIdList *list = cell->GetPointIds();

         if (c == VTK_TRIANGLE)
         {
             vtkIdType pt0 = list->GetId(0);
             vtkIdType pt1 = list->GetId(1);
             vtkIdType pt2 = list->GetId(2);
             for (int j = 0 ; j < nsteps-1 ; ++j)
             {
                 vtkIdType wedge[6];
                 wedge[0] = npts*j + pt0;
                 wedge[1] = npts*j + pt1;
                 wedge[2] = npts*j + pt2;
                 wedge[3] = npts*(j+1) + pt0;
                 wedge[4] = npts*(j+1) + pt1;
                 wedge[5] = npts*(j+1) + pt2;
                 if (j == nsteps-2 && overlap_ends)
                 {
                     wedge[3] = pt0;
                     wedge[4] = pt1;
                     wedge[5] = pt2;
                 }
                 ugrid->InsertNextCell(VTK_WEDGE, 6, wedge);
             }
         }
         else if(c == VTK_LINE)
         {
             vtkIdType pt0 = list->GetId(0);
             vtkIdType pt1 = list->GetId(1);
             vtkIdType npts_times_j = 0;
             for(int j = 0; j < nsteps - 1; ++j)
             {
                 vtkIdType quad[4];
                 if(j < nsteps-2)
                 {
                     quad[0] = npts_times_j + pt0;
                     quad[1] = npts_times_j + pt1;
                     quad[2] = npts_times_j + pt1 + npts;
                     quad[3] = npts_times_j + pt0 + npts;
                 }
                 else
                 {
                     quad[0] = npts_times_j + pt0;
                     quad[1] = npts_times_j + pt1;
                     quad[2] = pt1;
                     quad[3] = pt0;
                 }
                 ugrid->InsertNextCell(VTK_QUAD, 4, quad);
                 npts_times_j += npts;
             }
         }
         else if (c == VTK_VERTEX)
         {
             vtkIdType in_vert = list->GetId(0);
             vtkIdType npts_times_j = 0;
             for (int j = 0; j < nsteps - 1; ++j)
             {
                 vtkIdType out_vert[1];
                 if (j < nsteps-2)
                 {
                     out_vert[0] = npts_times_j + in_vert;
                 }
                 else
                 {
                     out_vert[0] = npts_times_j + in_vert;
                 }
                 ugrid->InsertNextCell(VTK_VERTEX, 1, out_vert);
                 npts_times_j += npts;
             }
         }
         else
         {
             vtkIdType pt0 = list->GetId(0);
             vtkIdType pt1 = list->GetId(1);
             vtkIdType pt2 = list->GetId(2);
             vtkIdType pt3 = list->GetId(3);
             if (c == VTK_PIXEL)
             {
                 pt2 = list->GetId(3);
                 pt3 = list->GetId(2);
             }
             for (int j = 0 ; j < nsteps-1 ; ++j)
             {
                 vtkIdType hex[8];
                 hex[0] = npts*j + pt0;
                 hex[1] = npts*j + pt1;
                 hex[2] = npts*j + pt2;
                 hex[3] = npts*j + pt3;
                 hex[4] = npts*(j+1) + pt0;
                 hex[5] = npts*(j+1) + pt1;
                 hex[6] = npts*(j+1) + pt2;
                 hex[7] = npts*(j+1) + pt3;
                 if (j == nsteps-2 && overlap_ends)
                 {
                     hex[4] = pt0;
                     hex[5] = pt1;
                     hex[6] = pt2;
                     hex[7] = pt3;
                 }
                 ugrid->InsertNextCell(VTK_HEXAHEDRON, 8, hex);
             }
         }
    }

    vtkCellData *incd   = in_ds->GetCellData();
    vtkCellData *outcd  = ugrid->GetCellData();
    outcd->CopyAllocate(incd, n_out_cells);
    for (vtkIdType i = 0 ; i < n_out_cells ; ++i)
    {
        outcd->CopyData(incd, i/(nsteps-1), i);
    }
    
    vtkPointData *inpd  = in_ds->GetPointData();
    vtkPointData *outpd = ugrid->GetPointData();
    outpd->CopyAllocate(inpd, n_out_pts);
    for (vtkIdType i = 0 ; i < n_out_pts ; ++i)
    {
        outpd->CopyData(inpd, i%npts, i);
    }
    
    //
    // Handle vector data.
    //
    vtkDataArray *vectors;
    vectors = ugrid->GetPointData()->GetVectors();
    if (vectors)
    {
        RevolveVectors(vectors, false, niter, npts, nsteps, 
                      stop_angle, start_angle, axis, mat, mt);
    }

    vectors = ugrid->GetCellData()->GetVectors();
    if (vectors)
    {
        RevolveVectors(vectors, true, nsteps-1, ncells, nsteps, 
                      stop_angle, start_angle, axis, mat, mt);
    }
 
    avtDataRepresentation *out_dr = new avtDataRepresentation(ugrid,
        in_dr->GetDomain(), in_dr->GetLabel());

    //
    // Clean up.
    //
    ugrid->Delete();
    mat->Delete();
    pts->Delete();

    return out_dr;
}


// ****************************************************************************
//  Method: avtRevolveFilter::UpdateDataObjectInfo
//
//  Purpose:
//      Tells the output that it is now higher in topological dimension.
//
//  Programmer: Hank Childs
//  Creation:   December 11, 2002
//
//  Modifications:
//    Kathleen Bonnell, Mon Apr 14 09:57:39 PDT 2003 
//    Set CanUseTransform to false.
//
//    Dave Bremer, Wed Dec 20 16:22:06 PST 2006
//    Set the coord system of the output to be XY (cylindrical meshes are
//    not cylindrical after you revolve them).
//
//    Hank Childs, Thu Aug 26 13:47:30 PDT 2010
//    Change extents names.
//
//    Brad Whitlock, Mon Apr  7 15:55:02 PDT 2014
//    Add filter metadata used in export.
//    Work partially supported by DOE Grant SC0007548.
//
// ****************************************************************************

void
avtRevolveFilter::UpdateDataObjectInfo(void)
{
    avtDataAttributes &inAtts      = GetInput()->GetInfo().GetAttributes();
    avtDataAttributes &outAtts     = GetOutput()->GetInfo().GetAttributes();
    avtDataValidity   &outValidity = GetOutput()->GetInfo().GetValidity();

    outAtts.SetTopologicalDimension(inAtts.GetTopologicalDimension()+1);
    if (inAtts.GetSpatialDimension() >= 2)
    {
        outAtts.SetSpatialDimension(3);
    }
    else
    {
        outAtts.SetSpatialDimension(inAtts.GetSpatialDimension()+1);
    }
    outValidity.InvalidateZones();
    outValidity.SetPointsWereTransformed(true);
    outValidity.InvalidateSpatialMetaData();

    //
    // This filter invalidates any transform matrix in the pipeline.
    //
    outAtts.SetCanUseTransform(false);

    //
    // The output will not be cylindrical.  Not cylindrical -> XY
    //
    outAtts.SetMeshCoordType(AVT_XY);

    //
    // Now revolve the extents.
    //
    double b[6];
    int olddim = inAtts.GetSpatialDimension();
    if (inAtts.GetOriginalSpatialExtents()->HasExtents())
    {
        inAtts.GetOriginalSpatialExtents()->CopyTo(b);
        RevolveExtents(b, olddim);
        outAtts.GetOriginalSpatialExtents()->Set(b);
    }

    if (inAtts.GetThisProcsOriginalSpatialExtents()->HasExtents())
    {
        inAtts.GetThisProcsOriginalSpatialExtents()->CopyTo(b);
        RevolveExtents(b, olddim);
        outAtts.GetThisProcsOriginalSpatialExtents()->Set(b);
    }

    if (inAtts.GetDesiredSpatialExtents()->HasExtents())
    {
        inAtts.GetDesiredSpatialExtents()->CopyTo(b);
        RevolveExtents(b, olddim);
        outAtts.GetDesiredSpatialExtents()->Set(b);
    }

    if (inAtts.GetActualSpatialExtents()->HasExtents())
    {
        inAtts.GetActualSpatialExtents()->CopyTo(b);
        RevolveExtents(b, olddim);
        outAtts.GetActualSpatialExtents()->Set(b);
    }

    if (inAtts.GetThisProcsActualSpatialExtents()->HasExtents())
    {
        inAtts.GetThisProcsActualSpatialExtents()->CopyTo(b);
        RevolveExtents(b, olddim);
        outAtts.GetThisProcsActualSpatialExtents()->Set(b);
    }

    outAtts.AddFilterMetaData("Revolve");
}


// ****************************************************************************
//  Method: avtRevolveFilter::RevolveExtents
//
//  Purpose:
//      Determines the extents of a dataset revolved around an axis.
//
//  Programmer: Hank Childs
//  Creation:   December 11, 2002
//
//  Modifications:
//
//    Hank Childs, Sun Mar 18 09:49:42 PDT 2007
//    Add support for ZR and RZ meshes.
//
// ****************************************************************************

void
avtRevolveFilter::RevolveExtents(double *dbounds, int spat_dim)
{
    int    nsteps      = atts.GetSteps();
    double start_angle = atts.GetStartAngle();
    double stop_angle  = atts.GetStopAngle();
    avtMeshCoordType mt = GetMeshType();
    double axis[3];
    GetAxis(mt, axis);

    vtkMatrix4x4 *mat = vtkMatrix4x4::New();
    double new_bounds[6];
    new_bounds[0] = FLT_MAX;
    new_bounds[1] = -FLT_MAX;
    new_bounds[2] = FLT_MAX;
    new_bounds[3] = -FLT_MAX;
    new_bounds[4] = FLT_MAX;
    new_bounds[5] = -FLT_MAX;

    //
    // Revolve the bounding box for each increment in theta.
    //
    for (int i = 0 ; i < nsteps ; ++i)
    {
        double angle = ((stop_angle-start_angle)*i)/(nsteps-1) + start_angle;
        GetRotationMatrix(angle, axis, mat, mt);
        int iters = 1;
        for (int k = 0 ; k < spat_dim ; ++k)
            iters *= 2;
        for (int j = 0 ; j < iters ; ++j)
        {
            double pt[4];
            pt[0] = (spat_dim >= 1 ? (j & 1 ? dbounds[1] : dbounds[0]) : 0.);
            pt[1] = (spat_dim >= 2 ? (j & 2 ? dbounds[3] : dbounds[2]) : 0.);
            pt[2] = (spat_dim >= 3 ? (j & 4 ? dbounds[5] : dbounds[4]) : 0.);
            pt[3] = 1.;

            double outpt[4];
            mat->MultiplyPoint(pt, outpt);
            new_bounds[0] = (outpt[0]<new_bounds[0] ? outpt[0] :new_bounds[0]);
            new_bounds[1] = (outpt[0]>new_bounds[1] ? outpt[0] :new_bounds[1]);
            new_bounds[2] = (outpt[1]<new_bounds[2] ? outpt[1] :new_bounds[2]);
            new_bounds[3] = (outpt[1]>new_bounds[3] ? outpt[1] :new_bounds[3]);
            new_bounds[4] = (outpt[2]<new_bounds[4] ? outpt[2] :new_bounds[4]);
            new_bounds[5] = (outpt[2]>new_bounds[5] ? outpt[2] :new_bounds[5]);
        }
    }

    mat->Delete();
    dbounds[0] = new_bounds[0];
    dbounds[1] = new_bounds[1];
    dbounds[2] = new_bounds[2];
    dbounds[3] = new_bounds[3];
    dbounds[4] = new_bounds[4];
    dbounds[5] = new_bounds[5];
}


// ****************************************************************************
//  Method: avtRevolveFilter::VerifyInput
//
//  Purpose:
//      Verifies that the input does not already have a topological dimension
//      of 3.
//
//  Programmer: Hank Childs
//  Creation:   December 11, 2002
//
// ****************************************************************************

void
avtRevolveFilter::VerifyInput(void)
{
    if  (GetInput()->GetInfo().GetAttributes().GetTopologicalDimension() >= 3)
    {
        EXCEPTION2(InvalidDimensionsException, "Revolve", "<=2D");
    }
}


// ****************************************************************************
//  Function: GetRotationMatrix
//
//  Purpose:
//      Given an angle and an axis, this creates a matrix that will rotate a
//      point around that axis.
//
//  Programmer: Hank Childs
//  Creation:   December 11, 2002
//
//  Modifications:
//
//    Hank Childs, Sun Mar 18 09:49:42 PDT 2007
//    Add support for RZ meshes.
//
// ****************************************************************************

static void
GetRotationMatrix(double angle, double axis[3], vtkMatrix4x4 *mat,
                  avtMeshCoordType mt)
{
    //
    // If we have an RZ mesh, then swap X and Z
    //
    vtkMatrix4x4 *rzCorrect = vtkMatrix4x4::New();
    rzCorrect->Identity();
    if (mt == AVT_RZ)
    {
        // SWAP X and Z
        rzCorrect->SetElement(0, 0, 0.);
        rzCorrect->SetElement(0, 2, 1.);
        rzCorrect->SetElement(2, 0, 1.);
        rzCorrect->SetElement(2, 2, 0.);
    }
    else if (mt == AVT_ZR)
    {
        // SWAP X and Y (to get to RZ), then swap X and Z
        // --> X->Y, Y->Z, Z->X
        rzCorrect->SetElement(0, 0, 0.);
        rzCorrect->SetElement(1, 1, 0.);
        rzCorrect->SetElement(2, 2, 0.);
        rzCorrect->SetElement(1, 0, 1.);
        rzCorrect->SetElement(2, 1, 1.);
        rzCorrect->SetElement(0, 2, 1.);
    }

    //
    // The game plan is to transform into a coordinate space that we are
    // familiar with, perform the rotation there, and then rotate back.
    //

    //
    // First rotate to the yz plane.  We will do this by rotating by theta
    // around the y-axis, where theta is arctan(axis[0] / axis[2]).
    //
    vtkMatrix4x4 *rot1 = vtkMatrix4x4::New();
    rot1->Identity();
    vtkMatrix4x4 *rot5 = vtkMatrix4x4::New();
    rot5->Identity();
    if (axis[0] != 0.)
    {
        double theta = atan2(axis[0], axis[2]);
        double cos_theta = cos(theta);
        double sin_theta = sin(theta);
        rot1->SetElement(0, 0, cos_theta);
        rot1->SetElement(2, 0, -sin_theta);
        rot1->SetElement(0, 2, sin_theta);
        rot1->SetElement(2, 2, cos_theta);

        //
        // Rot 5 will be a rotation around -theta.
        //
        double cos_minus_theta = cos_theta;
        double sin_minus_theta = -sin_theta;
        rot5->SetElement(0, 0, cos_minus_theta);
        rot5->SetElement(2, 0, -sin_minus_theta);
        rot5->SetElement(0, 2, sin_minus_theta);
        rot5->SetElement(2, 2, cos_minus_theta);
    }

    //
    // Now rotate around the x-axis until we get the to the z-axis.
    //
    vtkMatrix4x4 *rot2 = vtkMatrix4x4::New();
    rot2->Identity();
    vtkMatrix4x4 *rot4 = vtkMatrix4x4::New();
    rot4->Identity();
    if (axis[1] != 0.)
    {
        double theta = atan2(axis[1], sqrt(axis[0]*axis[0] + axis[2]*axis[2]));
        double cos_theta = cos(theta);
        double sin_theta = sin(theta);
        rot2->SetElement(1, 1, cos_theta);
        rot2->SetElement(2, 1, sin_theta);
        rot2->SetElement(1, 2, -sin_theta);
        rot2->SetElement(2, 2, cos_theta);

        //
        // Rot 4 will be a rotation around -theta.
        //
        double cos_minus_theta = cos_theta;
        double sin_minus_theta = -sin_theta;
        rot4->SetElement(1, 1, cos_minus_theta);
        rot4->SetElement(2, 1, sin_minus_theta);
        rot4->SetElement(1, 2, -sin_minus_theta);
        rot4->SetElement(2, 2, cos_minus_theta);
    }

    //
    // Now we can do the easy rotation around the z-axis.
    //
#if defined(_WIN32)
    double angle_rad = (angle / 360. * 2. * 3.14159);
#else
    double angle_rad = (angle / 360. * 2. * M_PI);
#endif
    vtkMatrix4x4 *rot3 = vtkMatrix4x4::New();
    rot3->Identity();
    double cos_angle = cos(angle_rad);
    double sin_angle = sin(angle_rad);
    rot3->SetElement(0, 0, cos_angle);
    rot3->SetElement(1, 0, sin_angle);
    rot3->SetElement(0, 1, -sin_angle);
    rot3->SetElement(1, 1, cos_angle);

    //
    // Now set up our matrix.
    //
    vtkMatrix4x4 *tmp  = vtkMatrix4x4::New();
    vtkMatrix4x4 *tmp2 = vtkMatrix4x4::New();
    vtkMatrix4x4 *tmp3 = vtkMatrix4x4::New();
    vtkMatrix4x4::Multiply4x4(rot5, rot4, tmp);
    vtkMatrix4x4::Multiply4x4(tmp, rot3, tmp2);
    vtkMatrix4x4::Multiply4x4(tmp2, rot2, tmp);
    vtkMatrix4x4::Multiply4x4(tmp, rot1, tmp3);
    vtkMatrix4x4::Multiply4x4(tmp3, rzCorrect, mat);

    tmp->Delete();
    tmp2->Delete();
    tmp3->Delete();
    rot1->Delete();
    rot2->Delete();
    rot3->Delete();
    rot4->Delete();
    rot5->Delete();
    rzCorrect->Delete();
}


